Chaotic Behavior in Shell Models and Shell Maps 
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We study the chaotic behavior of the "GOY" shell model by measuring the variation of the 
maximal Lyapunov exponent with the parameter e which determines the nature of the second 
invariant (the generalized "helicity" invariant). After a Hopf bifurcation, we observe a critical 
point at e c ~ 0.38704 above which the maximal Lyapunov exponent grows nearly linearly. For 
high values of e the evolution becomes regular again which can be explained by a simple analytic 
argument. A model with few shells shows two transitions. To simplify the model substantially we 
introduce a shell map which exhibits similar properties as the "GOY" model. 



Shell models of turbulence have been studied intensively by Kadanoff jl]-[| and many others (for a review see 
. It appears that a lot of properties of highly turbulent flows are nicely captured by those models and numerical 
computations of shell models are much more tractable than direct simulations of the Navier-Stokes equations. Of 
particular focus have been studies of intermittency effects where laminar quiescent periods are interrupted by strong 
. intermittency bursts of high energy dissipation. Most shell models follow completely deterministic dynamics but 
exhibits nevertheless strong intermittent behavior and this is believed to be caused by the internal chaotic dynamics. 
For the "GOY" shell model (|,||, the strength of this chaos is known to be related to the nature of the second invariant, 
i.e. the "helicity" invariant [jjj^]. It is the purpose of this paper to explore the nature of the chaotic behavior as 
the properties of the second invariant are changed. We do this by estimating the variation of the maximal Lyapunov 
^vq ' exponent with an external parameter e and find a transition in this exponent where it jumps up from zero to a 
small but finite value after which it grows linearly. Previous work by Biferale et al and Kadanoff et al || have 
concentrated on the instability of the "Kolmogorov fixed point" as the nature of the invariant is changed, measured 
by a parameter e to be introduced below. The first group observed a Hopf bifurcation of the fixed point at e = 0.3843 
turning into a torus at e = 0.3953 and finally a strange attractor at e = 0.398. This was refined by the other group 
who also studied the detailed eigenvalue spectrum for the unstable modes These studies however only addressed 
the instability of the stationary state and not the transition to chaotic evolution, which is our purpose here. 

Shell models are formed by various truncation techniques of the Navier-Stokes equations (J] . The most well-studied 
model is the "GOY" model of Gledzer-Ohkitani-Yamada |p]4l0|. This model yields corrections to the Kolmogorov 
theory Js) in good agreement with experiments []i~l| , [l2| . For the "GOY" shell model, wave-number space is divided into 
N separated shells each characterized by a wave- number k n — r" fco (f = 2), with n = 1, • • • N. Each shell is assigned 
. a complex amplitude u n describing the typical velocity gradient over a scale £ n — l/k n . By assuming interactions 
among nearest and next nearest neighbour shells and phase space volume conservation one arrives at the following 
evolution equations || 

^ ' d h r 

^~dt + n ' ^ = ifc "( a " M ™+l U «+2 + y W n-l M n+l + -J U n-lU n - 2 ) + /<W, (I) 

with boundary conditions b\ = b N = c\ = C2 = ajv_i = fljv = 0. / is an external, constant forcing, here on the forth 
mode. 

The coefficients of the non-linear terms must follow the relation a n + b n+ \ + c n+ 2 = in order to satisfy the 
conservation of energy, E = J2 n \ u n\ 2 , when / = v = 0. The constraints still leave a free parameter e so that one can 
set a n = l, b n+ i = — e, c n+ 2 = —(1 — e) 0. As observed by Kadanoff, one obtains the canonical value e = 1/2, if 
helicity conservation is also demanded Q . The set of N coupled ordinary differential equations can be numerically 
integrated by standard techniques. 

To compute the maximal Lyapunov exponent in the GOY model, we introduce the notation U = 
(i?e(ui), Im{ui), ■ • ■ , Re{uN), Im(u^) ) and Fi = dUi/dt and consider the linear variational equations 



* Electronic Address: julienk@sci.kun.nl 
^Electronic Address: okkels@nbi.dk 
'"Electronic Address: mhjensen@nbi.dk 



1 



, 2N 

^T=£A;-*j i = l,..,2N (2) 

i=i 

for the time evolution of an infinitesimal increment z = <5U, where 

A nj ee dFjdUj (3) 

is the Jacobian matrix of Eqs. ([lj). The solution for the tangent vector z can thus be formally written as zfa) — 
M(ti,t2) • with M = exp J^ 2 A(r)dT. A generic tangent vector z(£) is projected by the evolution along the 

eigenvector belonging to the maximum Lyapunov exponent, i.e. z(i) = |z(0)| e*- 1 ^ exp(Ai f) leading to 

Al = hm -lnM (4) 

t^oo t |z(0)| V ' 

where z(0) is the initial tangent vector. 

Practically, Eqs.[l], |^ are integrated simultaneously over a certain time St, starting with a normalized tangent vector 
in a random direction, z(0). The increment over time St in the length of the tangent vector is then Sz\ = |z(<5i)|/|z(0)|. 
Next, the tangent vector is normalized i{8t) = z(j£)/|z(<$t)| and this vector is used as a seed for a new integration 
over the time St i.e. propagated forward to t = 2St. Generalizing this argument we obtain the i'th increment 
Szi — \z(iSt)\/\z((i — l)St)\ and the maximal Lyapunov exponent is given by (where we now set A = Ai): 

M 



1 ^ ]n(6zj) 



A = lim 

Af->oo M ^— ' St 

t=l 

We have estimated the maximal Lyapunov exponent this way for the GOY model with standard parameters N — 
19, v = 1CP 6 , fco = 2~ 4 , / = (1 + i) * 0.005. The value of the time increment was varied from St — 0.001 all the way 
up to St =10, in time units set by the choice of parameters and we find that the results are independent of this value, 
as expected. Fig. 1 shows the obtained results, with a variation of e in the interval [0,2]. Note the transition to chaos 
at a critical value e c = 0.38704 [O. Biferale et al found a slightly larger value, e — 0.398, for the transition to 
the attractor at the same parameter values; we do not know the origin of this discrepancy. Magnifying around the 
transition point we observe a finite but small jump of A from up to 0.015. This makes the transition first order but 
we speculate that this jump might disappear in the N — > oo limit. After that the Lyapunov exponent grows more or 
less linearly with e and reaches a maximum value A ~ 0.54 at e = 0.92. Then it drops sharply down to zero around 
e ~ 1 after which it raises and drops again. The e = 1 is special because at this point the last term in the GOY model, 
which couples to the two previous shells, is zero. This means that for the last three shells, the equations reduce to: 

("^ + vk N-2) u N-2 = ik N - 2 {u* N _ x U* N - ^W^_ 3 li^_ 1 ) 

vk%_ x )u N -x = -ik N - 2 (u* N _ 2 u* N ) (6) 



y dt 

d 
dt 



It is seen that un decays exponentially to zero, and therefore we may neglect the right-hand-side of the second 
equation, which then leads to an exponential decay of itjv— 1< Just the same argument can be applied to the behavior 
of ujv_2 up to «4 where the forcing is applied. We thus expect that the solution essentially converges towards the 
trivial fixed point u n = ]14| ] and not to the "Kolmogorov fixed point" . This is in accordance with our numerical 
results, where we see a fast exponential convergence of the modes toward u n = 0. Fig. 2 shows the variation of the 
amplitudes \u n \ as a function of time. Besides the nice exponential decay, it appears that each mode is successively 
"triggered" for the decay. This is seen in Fig. 2: the n'th shell is triggered when the absolute value of the amplitude 
of the n + l'th shell is a approximately 5 • 10 4 and the amplitude of the n + 2'th shell has practically vanished. 
This means that the equation for the n'th shell becomes 4iU n — — i fc n _i(u*_ 1 u* +1 ) — vk^ u n , and just when the 
n + l'th amplitude reaches ~ 5 • 10 4 the right-hand-side will be dominated by the viscosity which then forces the n'th 
amplitude towards zero. 

Because the system is forced on the 4'th shell, the amplitude of this shell is prevented from going to zero. In fact, 
when the forcing is included the "trivial fixed point" is w* = (0, 0, 0, -Jgx, 0, 0, ). Fig. 3 shows the variation of the 
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real part of the maximal eigenvalue of the Jacobian as function of e evaluated in this fixed point. We observe that the 
fixed point is repelling for the standard parameter value e = 0.5 but becomes attractive for e = 1.0. This is also the 
place where the dynamics becomes non-chaotic as indicated by a convergence (although very slow) of the Lyapunov 
exponent towards zero. 

In order to understand the above results further we try to formulate the simplest possible model with the same 
"symmetries" as the GOY model, a shell map. This map is constructed from the same principles as the GOY model 
with the difference that the time interval between each update is a full time unit. We then obtain the following map 
in the real "velocity" variable v n (i) of shell n at integer time i: 

B C 

v n (i + l) =v n (i) + k n (A n v n+ i(i)v n+2 (i) + "y^n-i (i)v n +i (0 + v n -i(i)v n ~2(i)) - v\k 2 n v n (i) + /i<5 n ,i (7) 

with boundary conditions B\ = Bm = C\ = C2 = A/v-i = An = 0. Conservation of the energy E = ^2 n v n (i) 2 tells 
us that, for v\ = f\ = 0, the differential dE = 2J2 n v n (i)dv n (i) ~ 2^ n v„(i) • (v n (i + 1) — v n (i)) = leading to a 
similar relations between the coefficients as for the GOY model: A n = 1, B n = — ei, C n = — (1 — ei) Surprisingly 
enough, iterations of this model is stable even for very few shells. 

To simplify, we start out with a shell map with only TV = 5 shells and the following parameters: v = 10 , ko = 
2" 1 , / = 0.00005. It is immediate to study the eigenvalues of the Jacobian at the fixed point u n (i + 1) = u n (i). We 
locate the fixed point in the stable regime by iterating the map and refining the solution with Newton's method. This 
fixed point is a "Kolmogorov fixed point" with a scaling close to v„ ~ fc^ 1 ^ 3 , apart oscillations @. By varying the 
parameter t\ in small steps and then using Newton's method again one can obtain the fixed point where the map 
becomes unstable. We calculate the eigenvalues of the Jacobian at the fixed point and observe a Hopf-bifurcation (a 
pair of two complex conjugate eigenvalues escapes from the unit disc in the complex plane) at ei = 0.705081885. Here 
a limit cycle appears. From numerical iteration we observe that the limit cycle has a period of 217.822 (iterations) at 
ei = 0.70509. At this point, the phase of the eigenvalues equals 0.0288458 which is in agreement with the numerically 
found period since 2 • n/0. 0288458 = 217.820. Fig. 4a shows the limit cycle for the parameter value ei = 0.73 (note 
the somewhat unusual "return plot" where the n + l'th shell amplitude is plotted against the n'th shell amplitude. 
This gives nicer graphs that a usual return plot). As e% is increased further a series of standard period doubling 
bifurcations occur, Fig. 4b shows period four at t\ = 0.764. Finally, the the motion ends up on a strange attractor 
as shown in Fig. 4c for e\ = 0.77. 

To estimate the maximal Lyapunov exponent for the shell map a similar technique as for the GOY model is applied 
except that now we are dealing with a discrete map. The tangent vector z(i) is propagated according to 

N 

h{i + l) = ^A kj -Zj (8) 

where A k j is the Jacobian of the map @) 

A kj = dv k (i + l)/d Vj (i) (9) 

Again, practically we initiate by a unit tangent vector in a random direction z(0) and follow the expansion of the 
eigenvector. The only difference to the GOY model is that the time increment now is unity so Eq. (JsJ) applies with 
8t = l. 

Also at low values of t\ there is a transition and Fig. 5 show the variation of the maximal Lyapunov exponent 
as a function of t\. For e% ~ the model starts out to be chaotic but then the Lyapunov exponent drops and 
becomes negative around t\ ~ 0.15. The exponent approaches zero at ei ~ 0.705. We identify this transition with the 
appearance of the limit cycle: the distance between two initially close trajectories on the limit cycle will not converge 
nor diverge. Then the transition to chaos described above occurs at e\ c = 0.766. After that A grows sharply and 
drops back to zero at e% ~ 1. It is interesting that there seem to be windows where A goes to zero, which are like the 
windows found in the logistic map (see for instance [^6|). At e% = 1 we find again that the system converges towards 
the fixed point close to zero (^-,0,0,0,0) (and not the Kolmogorov fixed point). The eigenvalues of the Jacobian 

are in this case just 1 — vkf (i = 1, . . . 5) and thus the fixed point is stable. 

It is possible to apply the shell map also for more shells and one can go up to N — 11 without divergencies. Fig. 6 
shows the variation of the maximal Lyapunov exponent as a function of e^. Here the transition occurs at t\ c = 0.657 
For this number of shells one can identify an inertial range and below the transition, the structure functions follow 
the Kolmogorov predictions, although with the standard oscillations. Above e\ e the motion is intermittent and we 
again identify windows of stable evolution as for the N = 5 case. 
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FIG. 1. The variation of the maximal Lyapunov exponent as a function of e for the GOY model equations (jl|) with the 
parameters TV = 19, v — 10~ 6 , ko = 2~ 4 , / = (1 + i) * 0.005. Note the transition to chaos at e c = .38704 and the finite, but 
small jump at this point. The Lyapunov exponent increases nearly linearly but drops back almost to zero at e = 1 as explained 
in the text. The "bump" at e — 1.1 is due to the fact that the value of A has not converged. 

FIG. 2. A logarithmic plot of the long-term dynamics of the absolute value of the shell amplitudes |u n | for e = 1. The 
labelling refers to the shell number and one notes the drastic decay the of the higher shell amplitudes. 

FIG. 3. The value of the real part of the maximal eigenvalue of the Jacobian evaluated in the "trivial fixed point" 
— (0, 0, 0, -Jgt, 0, 0, ), as a function of e. Note, that the fixed point becomes stable around e = 1. 



FIG. 4. "Return plots" (v n (i),v n+1 (i)) of the shell map Eq. (0) with the parameters TV = 5, v = 10~ 4 ,fc = 2" 1 ,/ = 0.00005. 
In a) ei = 0.73 and we observe limit cycle generated by the Hopf bifurcation; in b) the second period doubling of the limit 
cycle is shown at ei = 0.764 and c) is the strange attractor for ei = 0.77. 



FIG. 5. The maximal Lyapunov exponent versus ei for the shell map in Eq. (0) for the parameters 
TV = 5, v = 10~ 4 ,fco = 2" 1 ,/ = 0.00005. Note the two transitions and that the model becomes regular again for e\ — 1. 



FIG. 6. The maximal Lyapunov exponent versus ei for the shell map in Eq. (|7|) for the parameters 
TV = 11, v — ko = 2~ 4 , / = 5 • 10~ 6 . Here, there is only one transition at ei c ~ 0.66 
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